import numpy as np
from scipy.signal import czt, find_peaks, filtfilt

# 控制干扰峰个数的 幅值判断
interference = 0.001


def get_czt(s, n, f1, f2, fs):
    """
    s:信号
    n:CZT点数
    f1:CZT变换后最低频率
    f2:CZT变换后最高频率
    fs:采样率
    return:
    返回CZT变换后的信号和频率
    z:CZT变换后归一化后的幅值 now_sig_czt
    fz:CZT变换后对应的频率 now_sig_freq
    """
    try:
        # 确保输入信号是数值数组
        s = np.asarray(s, dtype=np.float64)
        w = np.exp(-1j * 2 * np.pi * (f2 - f1) / (n * fs))
        a = np.exp(1j * 2 * np.pi * f1 / fs)
        z = czt(s, n, w, a)
        z = np.abs(z)
        fz = (np.arange(len(z)) * (f2 - f1) / len(z)) + f1
        z = z / np.max(z)
        fz = fz * 60
        return z, fz
    except Exception as e:
        print(f"Error in get_czt: {e}")
        return None, None


def ppg_process(channle_data, fs):
    """
    计算PPG,对信号进行处理
    """
    # 上面那个是对PPG信号做异常处理的，下面这个是做平滑的。
    abnormal_remove_ppg_channle = abnormal_remove(channle_data)
    abnormal_remove_ppg_channle = smooth_filt(abnormal_remove_ppg_channle, 6)

    filt_process_ppg_channle = filt_process(abnormal_remove_ppg_channle, fs)

    # 这个做动作干扰滤波的：
    motion_filt_process_ppg_channle = motion_filt_process(filt_process_ppg_channle)

    return motion_filt_process_ppg_channle


def acc_process(channle_data, fs):
    """
    计算ACC,对信号进行处理
    """
    # 上面那个是对ACC信号做异常处理的，下面这个是做平滑的。
    abnormal_remove_acc_channle = smooth_filt(channle_data, 6)

    filt_process_acc_channle = filt_process(abnormal_remove_acc_channle, fs)
    # 这个做动作干扰滤波的：
    motion_filt_process_acc_channle = motion_filt_process(filt_process_acc_channle)
    return motion_filt_process_acc_channle


def abnormal_remove(channle_data):
    """
    去除异常值
    """
    HF_a = [1, -3.590148110152597, 4.873443614170354, -2.963870062130241, 0.681973883426416]
    HF_b = [0.775707823645609, -3.094022833216837, 4.636644012395594, -3.094022833216837, 0.775707823645609]
    # 去均值
    channle_data = channle_data - np.mean(channle_data)
    # 一阶差分
    diff_channel_data = np.diff(channle_data)
    diff_channel_data = np.insert(diff_channel_data, 0, 0)
    # 找一阶差分波谷
    valley_index_list, _ = find_peaks(-diff_channel_data, distance=6)
    # 相邻两个波谷之间最大值位置找到被定义为波峰位置
    peak_index_list = []
    for temp_index in range(len(valley_index_list) - 1):
        peak_index = np.argmax(diff_channel_data[valley_index_list[temp_index]: valley_index_list[temp_index + 1]])
        peak_index += valley_index_list[temp_index]
        peak_index_list.append(peak_index)
    peak_index_list = np.array(peak_index_list)
    if len(peak_index_list) == 0:
        return channle_data
    # 波峰和前一个波谷的差的绝对值，以及波峰和后一个波谷的差的绝对值，放在一起排序，组成pks
    pks = np.abs(diff_channel_data[peak_index_list] - diff_channel_data[valley_index_list[0:-1]])
    pks = np.concatenate((pks, np.abs(
        diff_channel_data[peak_index_list] - diff_channel_data[valley_index_list[1:len(valley_index_list)]])))
    pks = np.sort(pks)
    if len(pks) == 0:
        return channle_data
    # 一阶差分去均值
    diff_channel_data -= np.mean(diff_channel_data)
    # 截止频率0.6，高通滤波
    diff_channel_data = filtfilt(HF_b, HF_a, diff_channel_data)
    # 通过差分信号重建ppg信号
    sig_list = [0]
    for temp_index in range(1, len(diff_channel_data)):
        sig_list.append(sig_list[-1] + diff_channel_data[temp_index])
    sig_list = np.array(sig_list)
    sig_list -= np.mean(sig_list)
    return sig_list


def smooth_filt(data, smooth_len):
    """
       去除异常值
       """
    result = []
    for index in range(len(data)):
        if index < len(data) - smooth_len:
            result.append(np.mean(data[index: index + smooth_len]))
    result.insert(0, result[0])
    result.insert(0, result[0])
    result.insert(0, result[0])
    result.append(result[-1])
    result.append(result[-1])
    result.append(result[-1])
    result = np.array(result)
    return np.array(result)


def filt_process(channle_data, FS):
    """
    信号滤波
    """
    HF_a = [1, -3.590148110152597, 4.873443614170354, -2.963870062130241, 0.681973883426416]
    HF_b = [0.775707823645609, -3.094022833216837, 4.636644012395594, -3.094022833216837, 0.775707823645609]
    LF_a = [1, -2.174805506915435, 2.357402921272809, -1.282544135294374, 0.306064117903952]
    LF_b = [0.036936953516246, 0.031977931425171, 0.061290122632184, 0.031977931425171, 0.036936953516246]
    # 首尾拼接一秒钟采样的数据
    temp_data = np.array(list(channle_data[FS:0:-1]) + list(channle_data) + list(channle_data[-2:-2 - FS:-1]))
    temp_data = temp_data - np.mean(temp_data)
    filtered_sig = filtfilt(HF_b, HF_a, temp_data)
    filtered_sig = filtfilt(LF_b, LF_a, filtered_sig)
    filtered_sig = filtered_sig[FS:-FS]
    return filtered_sig


def motion_filt_process(channle_data):
    """
    运动信号滤波
    """
    PF_b = [0.142537860520234, -0.215022543769845, 0.149319417312823, -0.215022543769845, 0.142537860520234]
    PF_a = [1, -2.67035772773509, 2.94351262911759, -1.63009314766181, 0.400438754415306]
    motion_filt_process_data = filtfilt(PF_b, PF_a, channle_data)
    return motion_filt_process_data


def read_numbers_from_file(file_path):
    numbers = []
    with open(file_path, 'r', encoding='utf-8-sig') as file:  # 使用 'utf-8-sig' 编码来自动处理BOM
        for line in file:
            line = line.strip()  # 去除换行符和空白字符
            if line:  # 确保行不为空
                try:
                    number = int(line)
                    numbers.append(number)
                except ValueError as e:
                    print(f"Error converting line to int: {line}. Error: {e}")
    return numbers


def determine_type(variable):
    if isinstance(variable, list):
        return "list"
    elif isinstance(variable, np.ndarray):
        return "numpy.ndarray"
    elif isinstance(variable, (int, float)):
        return "number"
    elif isinstance(variable, str):
        return "string"
    elif isinstance(variable, dict):
        return "dict"
    elif callable(variable):
        return "callable"
    else:
        return "unknown"


# def get_peak_of_wave(now_sig_czt):
# 	"""
# 	 找列表中的波峰
# 	:param now_sig_czt:
# 	:return: 返回波峰列表
# 	"""
# 	value_peak_index_list, _ = find_peaks(now_sig_czt, distance=6)
# 	return value_peak_index_list.tolist()
# 找列表中的波峰
def get_peak_of_wave(now_sig_czt, now_sig_freq):
    """
     找列表中的波峰,与波峰对应的幅值
    :param now_sig_czt:
    :return: value_peak_index_list:返回波峰数组;
    now_sig_freq:频率(心率值)数组
    now_sig_czt: 波峰对应的幅值数组
    """
    # 使用 find_peaks 找到信号中的波峰索引
    value_peak_index_list, _ = find_peaks(now_sig_czt, distance=6)

    # 获取对应的频率的位置
    frequency_peaks = now_sig_freq[value_peak_index_list]
    # 获取波形对应的幅值
    czt_values = now_sig_czt[value_peak_index_list]

    new_data = list(zip(value_peak_index_list, frequency_peaks, czt_values))
    out_data = []
    # #筛选波峰,只保留大于指定幅值(谱峰)的波峰
    for peak_data, fr_data, czt_data in new_data:
        if czt_data > 0.1:  # 指定幅值(czt_)
            out_data.append((peak_data, fr_data, czt_data))
    # 将 out_data 拆分成3个列表
    peak_list, frequency_peaks_list, czt_values_list = zip(*out_data) if out_data else ([], [], [])

    return peak_list, frequency_peaks_list, czt_values_list


# 找心率带波峰
def get_heart_rate_peak(now_sig_czt, now_sig_freq, now_polar_heartbeat):
    """
    找心率带波峰,对应的频率值,要求:最靠近心率峰与标签误差需小于10bpm,不满足这一条件时，这1s占比标记为0
    :param now_sig_czt:   z:CZT变换后归一化后的幅值 now_sig_czt
    :param now_sig_freq:   fz:CZT变换后对应的频率 now_sig_freq
    :param now_polar_heartbeat: 金标心率
    :return:返回对应最接近金标心率的波峰的 对应的幅值(czt_value),否则返回0
    """

    # 找列表中的波峰
    value_peak_index_list, _ = find_peaks(now_sig_czt, distance=6)
    # 频谱中,心率带心率的位置
    polar_heart_rate_nearest_index = \
        np.where(np.abs(now_sig_freq - now_polar_heartbeat) == np.min(np.abs(now_sig_freq - now_polar_heartbeat)))[0][0]
    # 寻找离心率峰最近的频谱峰的位置
    polar_hr_peaks_diff = np.abs(np.array(value_peak_index_list) - polar_heart_rate_nearest_index)
    if len(polar_hr_peaks_diff) == 0:
        nearest_polar_czt_index = 0
        return nearest_polar_czt_index
    else:
        nearest_polar_peak_index = np.where(polar_hr_peaks_diff == np.min(polar_hr_peaks_diff))[0][0]
        nearest_polar_czt_index = value_peak_index_list[nearest_polar_peak_index]
        # 获取对应的频谱峰的位置
        frequency_value = now_sig_freq[nearest_polar_czt_index]
        # 获取对应的幅值的位置
        czt_value = now_sig_czt[nearest_polar_czt_index]

        if abs(frequency_value - now_polar_heartbeat) < 10:
            return czt_value
        else:
            return 0


# 最大峰能量占比
def max_peak_energy_ratio(czt_values):
    """
    最大峰能量占比
    :param czt_values:幅值列表
    :return: 返回占比值
    """
    # 获取幅值最大值
    czt_values_max = np.max(czt_values)
    print("幅值最大值:" + str(czt_values_max))
    czt_values_sum = sum(czt_values)
    print("幅值合:" + str(czt_values_sum))
    return czt_values_max / (sum(czt_values))


# 计算干扰峰个数
def interfering_peaks_count(now_sig_czt, now_sig_freq):
    """
    计算干扰峰个数,要求:波峰大于0.5则为干扰峰
    :param peak_list: 波峰列表
    :return:返回干扰峰个数
    """

    # 使用 find_peaks 找到信号中的波峰索引
    value_peak_index_list, _ = find_peaks(now_sig_czt, distance=6)

    # 获取波形对应的幅值
    czt_values = now_sig_czt[value_peak_index_list]

    # 获取对应的频率的位置
    frequency_peaks = now_sig_freq[value_peak_index_list]

    # 计算peak_list中大于0.5的个数
    # 把peak_list中大于0.5的数添加进一个新列表
    peak = []
    for i in czt_values:
        if i > interference:
            peak.append(i)

    return len(peak)


# 心率峰能量占比
def energy_ratio(czt_values_list, czt_values):
    """
    计算能量占比
    :param frequency_peaks:幅值数组
    :param czt_values: 金标所在索引位置对应的频率值 金标对应的幅值
    :return:返回心率峰能量占比
    """
    # 计算能量占比
    # print("频心率峰能量合:",str(sum(frequency_peaks)))
    energy_ratio = czt_values / sum(czt_values_list)

    return energy_ratio


def determine_type(variable):
    """
    确认类型,用于被c调用时确认
    :param variable:
    :return:字符串
    """
    if isinstance(variable, list):
        return "list"
    elif isinstance(variable, np.ndarray):
        return "numpy.ndarray"
    elif isinstance(variable, (int, float)):
        return "number"
    elif isinstance(variable, str):
        return "string"
    elif isinstance(variable, dict):
        return "dict"
    elif callable(variable):
        return "callable"
    else:
        return "unknown"


# 计算AC/DC能量比值
def get_ac_dc_energy_ratio(window_data, Gain_list):
    """
    计算AC/DC能量比值
    :param window_data: 未处理前的数据,10s窗口的数据
    :param Gain_list: 增益列表
    :return:返回AC/DC能量比值
    """

    # 设置增益
    gain_value = np.array([10, 25, 50, 75, 100, 250, 500, 750, 1000]) / 5

    # 处理增益值
    gain_data_list = gain_value[np.array(Gain_list)]

    # 计算直流能量（DC energy）
    dc_energy = (np.array(window_data) - 2 ** 23) / gain_data_list

    # 计算交流能量（AC energy）
    x = np.array(window_data)
    # 计算信号的均值
    mean_x = np.mean(x)
    # 去掉均值
    x_prime = x - mean_x
    # 使用高通滤波
    ac_energy = high_pass_filtering(x_prime, 25)

    # ac/dc 能量比值列表
    ac_dc_energy_list = ac_energy / dc_energy

    # 计算 25% 和 75% 分位数
    q25 = np.percentile(ac_dc_energy_list, 25)
    q75 = np.percentile(ac_dc_energy_list, 75)
    return (q25 + q75) / 2


def high_pass_filtering(channle_data, FS):
    """
    信号进行高通滤波
    """
    HF_a = [1, -3.590148110152597, 4.873443614170354, -2.963870062130241, 0.681973883426416]
    HF_b = [0.775707823645609, -3.094022833216837, 4.636644012395594, -3.094022833216837, 0.775707823645609]

    # 首尾拼接一秒钟采样的数据
    temp_data = np.array(list(channle_data[FS:0:-1]) + list(channle_data) + list(channle_data[-2:-2 - FS:-1]))
    temp_data = temp_data - np.mean(temp_data)
    # 对信号进行双向滤波
    filtered_sig = filtfilt(HF_b, HF_a, temp_data)
    filtered_sig = filtered_sig[FS:-FS]
    return filtered_sig


# 获取频域相关数据
def get_cztdata_process(window_data, fs, Gain_list=None, now_polar_heartbeat=None):
    """
           获取频域相关数据
    :param window_data: 窗口数据(ppg光路原始数据)
    :param fs: 采样率
    :param Gain_list: 增益列表 用于计算ac/dc能量比值 默认 = None
    :param now_polar_heartbeat: 金标心率 默认 = None
    :return:
           now_sig_czt:幅值
           now_sig_freq:频率(心率值)数组
           max_peak_of_wave:最大峰能量占比
           some_peak_of_wave:干扰峰个数
           energy_ratio_of_wave:心率峰能量占比
           peak_list:波峰数组
           frequency_peaks_list:波峰对应的频率数组
           czt_values_list:波峰对应幅值数组
           ac_dc_energy:ac/dc能量占比
    """
    # 默认心率峰能量占比返回0
    energy_ratio_of_wave = 0
    max_peak_of_wave = 0
    some_peak_of_wave = 0
    ac_dc_energy = 0
    # 1.预处理滤波等
    data = ppg_process(window_data, fs)
    # 2.先进行ctz转换
    now_sig_czt, now_sig_freq = get_czt(data, 256, 0.6, 4, 25)
    # 3.再进行处理
    # 找指定的大于指定幅值 波峰列表
    peak_list, frequency_peaks_list, czt_values_list = get_peak_of_wave(now_sig_czt, now_sig_freq)

    # 最大峰能量占比
    max_peak_of_wave = max_peak_energy_ratio(czt_values_list)
    # 干扰峰个数
    some_peak_of_wave = interfering_peaks_count(now_sig_czt, now_sig_freq)
    if Gain_list is not None:
        # 获取AC/DC能量比值
        ac_dc_energy = abs(get_ac_dc_energy_ratio(window_data, Gain_list))

    if now_polar_heartbeat is not None:
        if now_polar_heartbeat > 0:
            hr = get_heart_rate_peak(now_sig_czt, now_sig_freq, now_polar_heartbeat)
            # print("心率波峰位置=" + str(hr))
            # 心率峰能量占比
            energy_ratio_of_wave = energy_ratio(czt_values_list, hr)
    # print("心率峰能量占比=" + str(energy_ratio_of_wave))
    print("最大峰能量占比=" + str(max_peak_of_wave))
    print("干扰峰个数=" + str(some_peak_of_wave))
    print("心率峰能量占比=" + str(energy_ratio_of_wave))
    print("ac/dc能量比值=" + str(ac_dc_energy))
    return now_sig_czt, now_sig_freq, max_peak_of_wave, some_peak_of_wave, energy_ratio_of_wave, peak_list, frequency_peaks_list, czt_values_list, ac_dc_energy

